// Copyright (c) 2009-2010, Jeremy Brewer
// All rights reserved.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are met:
//
// * Redistributions of source code must retain the above copyright notice,
// this list of conditions and the following disclaimer.
// * Redistributions in binary form must reproduce the above copyright
// notice, this list of conditions and the following disclaimer in the
// documentation and/or other materials provided with the distribution.
//
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
// POSSIBILITY OF SUCH DAMAGE.

package galaxie500.numeric;

import java.util.Arrays;

import galaxie500.exceptions.NotImplementedError;
import galaxie500.exceptions.ValueError;
import galaxie500.util.MathUtil;

/**
 * Class implementing monotonic cubic spline interpolation
 * 
 * This class performs interpolation using cubic splines such that the input
 * function values being interpolated are guaranteed to be the extrema of the
 * interpolated function. In other words, the splines between the points are
 * monotonically increasing or decreasing.
 * 
 * The interpolation used is also robust to sensitive adjustments of the
 * function values in cases where extremely sharp variations (e.g. step
 * functions) are present.
 * 
 * The algorithm used is taken from an astrophysics paper:
 * Steffen, M., Astro. & Astrophysics, 239, 443-450, 1990
 * 
 * @author Jeremy Brewer
 */
public class Interpolator implements Function {

  /** Array of x values, copied by reference from the array passed to the
   * constructor.
   */
  private double[] xValues = null;

  /** Array of cubic polynomials, same size as xValues. The polynomials are
   * evaluated in x - xValues[i] where i is such that
   * xValues[i] <= x < xValues[i + 1].
   */
  private Polynomial[] polynomials = null;

  public Interpolator(double[] x, double[] y) {
    assert (x.length > 0);
    assert (x.length == y.length);
    this.xValues = x;
    polynomials = new Polynomial[x.length];
    computeCoefficients(y);
  }

  public static double integrate(double[] x, double[] y, double a, double b) {
    Interpolator i = new Interpolator(x, y);
    return i.integrate(a, b);
  }

  private void computeCoefficients(double[] y) {
    // TODO: Deal with corner cases here    
    if (xValues.length < 4) {
      throw new NotImplementedError(
          "Corner case for small (< 4) number of points not implemented");
    }
    
    double[] dydx = new double[xValues.length];
    double[] p = new double[xValues.length];
    double[] h = new double[xValues.length - 1];
    double[] s = new double[xValues.length - 1];
    
    // Calculate the derivatives on the grid points.
    h[0] = xValues[1] - xValues[0];
    s[0] = MathUtil.safeDivide(y[1] - y[0], h[0]);

    for (int i = 1; i < xValues.length - 1; i++) {
      h[i] = xValues[i + 1] - xValues[i];
      s[i] = MathUtil.safeDivide(y[i + 1] - y[i], h[i]);
      p[i] = MathUtil.safeDivide(s[i - 1] * h[i] + s[i] * h[i - 1],
                                 h[i - 1] + h[i]);
      dydx[i] = (Math.signum(s[i - 1]) + Math.signum(s[i])) *
                MathUtil.min3(Math.abs(s[i - 1]),
                              Math.abs(s[i]),
                              0.5 * Math.abs(p[i]));
    }

    // Makes expressions below much easier to read.
    int len = xValues.length;

    double tmp = MathUtil.safeDivide(h[0], h[0] + h[1]);
    p[0] = s[0] * (1.0 + tmp) - s[1] * tmp;
    double tmp2 = MathUtil.safeDivide(h[len - 2], h[len - 2] + h[len - 3]);
    p[len - 1] = s[len - 2] * (1.0 + tmp2) - s[len - 3] * tmp2;

    dydx[0] = (Math.signum(s[0]) + Math.signum(p[0])) *
              Math.min(Math.abs(s[0]), 0.5 * Math.abs(p[0]));
    dydx[len - 1] = (Math.signum(s[len - 2]) + Math.signum(p[len - 1])) *
                    Math.min(Math.abs(s[len - 2]), 0.5 * Math.abs(p[len - 1]));

    // Calculate the coefficients of the interpolating function.
    for (int i = 0; i < xValues.length - 1; i++) {
      double[] coef = new double[4];      
      coef[3] = MathUtil.safeDivide(dydx[i] + dydx[i + 1] - 2.0 * s[i],
                                    MathUtil.square(h[i]));
      coef[2] = MathUtil.safeDivide(3.0 * s[i] - 2.0 * dydx[i] - dydx[i + 1],
                                    h[i]);
      coef[1] = dydx[i];
      coef[0] = y[i];
      polynomials[i] = new Polynomial(coef);
    }

    // These extra coefficients are only used for the case where exact value of
    // the upper boundary is requested.
    polynomials[len - 1] =
        new Polynomial(new double[] {y[len -1], dydx[len - 1], 0.0, 0.0});
  }

  /**
   * Returns whether x can be interpolated.
   * 
   * @param x The value to test
   * @return Whether x is contained within the valid interpolation interval
   */
  public boolean contains(double x) {
    return x >= xValues[0] && x <= xValues[xValues.length - 1];
  }
  
  /**
   * Returns i such that xValues[i] <= x < xValues[i + 1].  If x lies outside
   * of xValues, returns -1.  If x == xValues[xValues.length - 1], returns
   * xValues.length - 1.
   * 
   * @param x The point to search for
   * @return The index as described above
   */
  private int bisect(double x) {
    if (!contains(x)) return -1;
    int i = Arrays.binarySearch(xValues, x);
    if (i < 0) return -(i + 2);
    return i;
  }
  
  /**
   * Returns the interpolated value at x.  If x lies outside of the valid
   * interpolation region, throws a ValueError.
   * 
   * @see galaxie500.numeric.Function#valueAt(double)
   * @throws ValueError If x lies outside of the valid x range
   */
  @Override
  public double valueAt(double x) {
    // Out of bounds check.
    int i = bisect(x);
    if (i < 0) {
      throw new ValueError(invalidPointError(x));
    }
    return polynomials[i].valueAt(x - xValues[i]);
  }

  /**
   * Returns an error message for a point lying outisde of the valid
   * interpolation region.
   * 
   * @param x The bad point
   * @return A string error message displaying the point and valid region
   */
  private String invalidPointError(double x) {
    return String.format(
        "Value %.14e is outside of array values [%.14e, %.14e]",
         x, xValues[0], xValues[xValues.length - 1]);
  }
  
  /**
   * Integrates the interpolated function over [a, b].
   * 
   * @throws ValueError If a or b lies outside the valid interpolation region
   * @param a The starting point of integration
   * @param b The end point of integration
   * @return The integral from a to b
   */
  public double integrate(double a, double b) {
    int i = bisect(a);
    if (i < 0) {
      throw new ValueError(invalidPointError(a));
    }

    int j = bisect(b);
    if (j < 0) {
      throw new ValueError(invalidPointError(b));
    }

    double integral = 0.0;

    // Sum up interior points.  Recall that the polynomial is in x - xValues[i].
    for (int k = i + 1; k < j; k++) {
      integral += polynomials[k].integrate(0.0, xValues[k + 1] - xValues[k]);
    }
    
    // Sum up leftover for end points.
    integral += polynomials[i].integrate(a - xValues[i],
                                         xValues[i + 1] - xValues[i]);
    integral += polynomials[j].integrate(0.0, b - xValues[j]);

    return integral;
  }
}
